//***************************************
//*         SIMPLEX METHOD for {LP       *
//*            BIG M  APPROACH          *
//*  A.D.BELEGUNDU & T.R.CHANDRUPATLA   *
//***************************************

dim pnp, nv, nl, ng, ne, nr, nc, li
dim a(), b(), nbs()

private sub simplexlp()
     bigm = 0
     for {j = 0; j < nv; j++) {
        bigm = bigm + pnp * Math.abs(a[nr][j])
     }
  //-----  slack, artificial and surplus variables
     if (nl > 0 ) {
        for {i = 0; i < nl; i++) {
           a[i][nv + i] = 1
           nbs[i] = nv + i
        }
     }
     if (ng > 0 ) {
        for {i = 0; i < ng; i++) {
           a[nl + i][nv + nl + i] = -1
           a[nl + i][nv + nl + ng + i] = 1
           nbs[nl + i] = nv + nl + ng + i
        }
     }
     if (ne > 0 ) {
        for {i = 0; i < ne; i++) {
           a[nl + ng + i][nv + nl + 2 * ng + i] = 1
           nbs[nl + ng + i] = nv + nl + 2 * ng + i
        }
     }
        nge = ng + ne
     if (nge > 0 ) {
        for {i = 0; i < nge; i++) {
           a[nr][nv + nl + ng + i] = bigm
        }
     }
     iter = 0
     printtableau()
  //----- removal of artificial variables from the objective row
     if (nge > 0 ) {
        for {i = 0; i < nge; i++) {
           c = a[nr][nv + nl + ng + i]
           for {j = 0; j < nc; j++) {
              a[nr][j] = a[nr][j] - c * a[nl + i][j]
           }
           b[nr] = b[nr] - c * b[nl + i]
        }
        printtableau()
     }
  //----- simplex routine
  //----- check criterion
     do {
     ict = 0
        for {jj = 0; jj < nc; jj++){
           c = a[nr][jj]
           if (c < 0 ) {
              ict = ict + 1
              if (ict = 1 ) {
                 amin = c
				 icv = jj
              }else{
                 if (amin > c ) {
                    amin = c
					icv = jj
                 }
              }
           }
        }
        if (ict == 0 ) {break}
  //----- incoming variable is icv
  //pivot row ipt ==> variable leaving basis
        ict = 0
        for {i = 0; i < nr - 1; i++) {
           if (a[i][icv] > 1e-30 * (1 + b[i]) ) {
              ict = ict + 1
              c1 = b[i] / a[i][icv]
              if (ict == 1 ) {
                 c = c1
				 ipt = i
              }else{
                 if (c1 < c ) {
                    c = c1
					ipt = i
                 }
              }
           }
        }
        nbs[ipt] = icv
        if (ict == 0 ) {
           Dummy  =  Dummy + "objective fn. is unbounded..."
           break
        }
  //----- basis change operation
        c1 = 1 / a[ipt][icv]
		b[ipt] = c1 * b[ipt]
        for {j = 0; j < nc; j++) {
           a[ipt][j] = c1 * a[ipt][j]
        }
        for {i = 0; i < nr; i++) {
           if (i <> ipt ) {
              c2 = a[i][icv]
              for {j = 1 to nc
                 a[i][j] = a[i][j] - c2 * a[ipt][j]
              next j
              b[i] = b[i] - c2 * b[ipt]
           }
        }
        iter = iter + 1
        printtableau()
     } while (ic < 1)
end sub
function printtableau()
        //----- print tableau
        li = li + 2
     if (iter = 0 ) {
        Tableu = Tableu + "initial/adjusted tableau" "'\n'
     }else{
        Tableu  = Tableu + "iter = " + iter + " incvar = " + icv + " pivrow = " + ipt + '\n'
     }
        for {i = 0; i < nr; i++) {
           if (i < nr ) {
              Tableu = Tableu + nbs[i]
           }else{
              Tableu = Tableu +  "obj "
           }
           for (j = 0; j < nc; j++) {
              Tableu = Tableu +  a[i][j]
           }
           Tableu = Tableu +  b[i]
        }
end sub
private sub output()
  //----- output on sheet 1

     Dummy = Dummy + "solution (see tableau at the bottom)" + '\n'
     Dummy = Dummy +  "variable  value" + '\n'
     for {i = 0; i < nr - 1; i++) {
        Dummy = Dummy +  nbs[i] + '\n'
        Dummy = Dummy +  b[i] + '\n'
     }
     Dummy = Dummy +  "maximum function value = " + '\n'
     Dummy = Dummy +  b[nr] + '\n'
end sub


